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RAPID CALCULATION OF THEORETICAL CMB ANGULAR POWER SPECTRA 



(N 
O 
O 
(N 

< 



> 
m 

o 

O 

^: 

Or 

6 ■ 

£: 



X 



Manoj Kaplinghat, Lloyd Knox and Constantinos Skordis 

Department of Physics, One Shields Avenue 
University of California, Davis, California 95616, USA 
Submitted to ApJ 

ABSTRACT 

We have developed a fast method for predicting the angular power spectrum, Ci, of the cosmic mi- 
crowave background given cosmological parameters and a primordial power spectrum of perturbations. 
After pre-computing the radiation temperature and gravitational potential transfer functions over a 
small sub-space of the total model parameter space, the rest of the model space (six or more cosmo- 
logical parameters and arbitrarily many primordial power spectrum parameters) is reached via rapid 
analytic and semi-analytic approximations which are highly accurate on all angular scales for which lin- 
ear perturbation theory applies. A single power spectrum can be calculated in ~ 1 second on a desktop 
computer. We discuss applications to cosmological parameter estimation. 

Subject headings: cosmology: theory - cosmic microwave background 



1. INTRODUCTION 

The anisotropy of the Cosmic Microwave Background 
(CMB) is proving to be a powerful cosmological probe. 
Measurements of its angular power spectrum can be used 
to tell us about the baryon density, dark energy density, 
the nature of the dark matter, the age of the universe and 
the primordial spectrum of perturbations generated in the 
inflationary era (Pryke et al. 2001; Netterfield et al. 2001; 
Lee et al. 2001; Wang et al. 2001; Knox et al. 2001). 

A persistent challenge to the analysis is the large num- 
ber of model angular power spectra (Cj) that must be 
calculated in order to understand the constraints the data 
place on parameter spaces with seven to ten or even higher 
dimensions. Here we present a fast, yet accurate, method 
for computing the Ci for a given model. 

These model angular power spectra are the expectation 
values of the variance of spherical harmonic coefficients, 
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The Cj depend on 



the density of dark matter, the fraction of this which is hot 
dark matter, the density of baryonic matter, the redshift of 
reionization of the intergalactic medium, the dark energy 
density, the dark energy pressure and the mean spatial 
curvature. These cosmological parameters influence the 
evolution of perturbations in the photon temperature. The 
Ci 's also depend on the statistical properties of the initial 
perturbations, possibly produced in an epoch of inflation. 
These initial conditions are described with the primordial 
gravitational potential power spectrum, P(k). 

The Ci can be calculated highly accurately because of 
the applicability of linear perturbation theory. Indeed, this 
is one of the reasons the CMB is such a powerful cosmo- 
logical probe. One need only solve the linearized Einstein 
and relevant Boltzmann equations, which can be cast as a 
set of coupled ordinary differential equations. Early codes 
(e.g. Bond & Efstathiou (1984)) directly solved the whole 
hierarchy, up to some limiting multipole moment, of the 
photon temperature perturbation and could take tens of 
hours to calculate C; for a single model. 

Hu & Sugiyama (1995) introduced a semi-analytic ap- 
proach which was much faster than the "whole hierarchy" 
Boltzmann codes of the day and with an accuracy around 
10%. That accuracy could be improved, but only at the 



expense of much slower performance. 

The line-of-sight integration method for solving the lin- 
earized Einstein and Boltzmann equations (Seljak & Zal- 
darriaga 1996) greatly reduced the time required for cal- 
culation of accurate theoretical power spectra by bypass- 
ing the need to solve the whole hierarchy. Publicly avail- 
able codes based on this method, mostly CMBfast (Sel- 
jak & Zaldarriaga 1996), have been the workhorses of all 
parameter-determination efforts to date. Despite its great 
speed, these analysis efforts have typically required months 
of running CMBfast. 

Tegmark & Zaldarriaga (2000) introduced a high-€ / 
\ow-£ split in the calculation of Ci to exploit analytic ap- 
proximations valid at high I and insensitivity to certain 
parameters at low I. With this split they were able to cal- 
culate a 7-dimensional grid of C; 's with many fewer calls of 
CMBfast than would have been required for a brute-force 
calculation. By analytically correcting these for reioniza- 
tion effects and scaling the tensor and scalar power spectra 
with separate amplitudes they covered a 10-dimensional 
parameter space. 

We use a similar high-^ / low-£ split to exploit the same 
analytic approximations as Tegmark & Zaldarriaga (2000). 
Our methods though have several advantages including 
reduced pre-compute time (~ 1 month reduced to ~ 40 
hours), and ability to handle large numbers of primordial 
power spectrum parameters. We have performed exten- 
sive accuracy tests showing that our errors are smaller 
than cosmic variance errors for I < 1000. We achieve 
these advantages by storing the Fourier- and Legendre- 
transformed temperature perturbation transfer function 
(rather than C;), using more efficient choices for grid pa- 
rameters, and further use of (highly accurate) analytic ap- 
proximations. We also have an option where all the \ow-£ 
effects are calculated with semi-analytic approximations, 
greatly reducing pre-compute times and storage require- 
ments even further, as well as allowing for greater ease in 
incorporating new physical effects. Our software package 
is called the Davis Anisotropy Shortcut (DASh) 1 . 
i 

DASh can be downloaded from 

http: / /www. physics. ucdavis.edu/Cosmology/dash/ 
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DASh incorporates many approximations that have been 
presented elsewhere in the literature, and also some new 
ones. We present an approximate scaling for the polariza- 
tion \ow-£ reionization feature, an improved approximation 
for the reionization damping factor and a generalization of 
the angular-diameter distance scaling which makes it ac- 
curate for all angular scales even in the presence of non- 
zero curvature. We also present improved semi-analytical 
approximations for the calculation of low-£ temperature 
spectra. These improvements enable a semi-analytical cal- 
culation of low-£ spectra to an accuracy better than 2% on 
average over a wide range of cosmological parameters (in- 
cluding curvature). 

Although recent work has shown the exploration of these 
large model spaces to be possible without DASh, our method 
greatly reduces the required computer resources. As such 
it will allow for extension to more parameters including 
those needed to describe isocurvature components, or the 
dark energy pressure. A particulary straightforward ex- 
tension would be to the number of parameters used to 
describe P(k), beyond the usual two needed for the power- 
law description. A faster method also makes it possible to 
redo calculations to check for sources of systematic error. 

A preliminary version of DASh has already been used 
for parameter estimation from CMB data (Knox et al. 
2001). There we combined DASh with the Monte Carlo 
Markov Chain (MCMC) approach to Bayesian inference 
described in Christensen et al. (2001). The MCMC ap- 
proach requires many fewer likelihood evaluations than 
a direct grid-based approach even for applications with 
only a handful of parameters, and generally becomes even 
more advantageous as the dimensionality increases further 
(Gilks et al. 1996). Others have used MCMC for cosmolog- 
ical problems (e.g. Verde & Spergel 2002) and we expect 
the technique to become widely used in cosmology. 

CMB anisotropies are conveniently broken up into two 
different types: those which are simply projections of fea- 
tures on (or near) the last-scattering surface (early anisotropy) 
and those that are generated much more recently (late 
anisotropy). After reviewing some notation in Section 2 
we discuss our computation of early anisotropies in Sec- 
tion 3. In Section 4 we describe the two different ways we 
compute late-time effects due to gravitational potential 
decay and reionization of the inter-galactic medium. In 
Section 5 we quantify the level of accuracy by comparing 
6,823 models as calculated with DASh to those calculated 
with CMBfast. In section 6 we describe our calculation of 
polarization power spectra and the contribution from ten- 
sor perturbations. In Section 7 we consider extensions, for 
example to including lensing effects, and finally in Section 
7 we conclude. 

2. NOTATION 

Before discussing the method we quickly review some 
notation. The temperature observed in direction 7 ob- 
served from any point in space, x, can be written as 

T(x, 7 )=T + A(x,7). (1) 
For anisotropy sourced by scalar metric perturbations the 
Fourier transformed temperature perturbation is azimuthally 
symmetric and can be expanded in Legendre polynomials 
as 

A(k, 7) = £ (21 + 1) (-ij'AKk^Oi) (2) 



where n = k ■ 7. The multipole moments of the Fourier- 
transformed temperature perturbation can be written as 
A/(k) = *j(k)A;(fc) where k = kk and *i(k) is the per- 
turbation in the gravitational potential (Ma & Bcrtschinger 
1995) at some very early time when all relevant perturba- 
tion wavelengths are larger than the horizon. Note that 
when we write Aj (k) with a scalar rather than vector argu- 
ment (as we do throughout), we are using it as a transfer 
function. 

If we solve for A;(fc) assuming adiabatic initial condi- 
tions with \ti(k) = 1, then if we assume the perturbations 
are statistically isotropic and homogeneous we can calcu- 
late Ci for any arbitrary initial potential power spectrum 
P(k) as 

/•OO 

d = (4tt) 2 / dkk 2 A 2 (k)P(k) (3) 
Jo 

where C; is defined by 

(a; m (x)a*, TO ,(x)) =CiS w S mm/ (4) 

and 

a im (x) = J d7Y; m (7)A(x,7). (5) 

We often express densities in units of the critical density 
for h = 1 where H = lOOh kmsec -1 Mpc -1 and the criti- 
cal density is p c = 3Hq / (8ttG) . Following convention, we 
refer to densities in these units with the symbol lo. The 
baryon density is ujb, the dark matter density is u>d, the 
matter density is u> m = u>t, + Ud, and the dark energy den- 
sity is cj x . Note that u>i = Vliti 2 . These symbols all refer 
to present day densities. We define a curvature "density" 
as lok = Qi{h 2 = (1 — f2tot)fr 2 where K = — 1,+1,0 cor- 
responds to an open, closed or flat universe respectively. 
With this definition the Friedmann equation at the present 
time becomes h 2 = We assume that a fraction, 

fh, of the dark matter is hot and that the rest is cold. 
We further assume that the dark energy is a cosmologi- 
cal constant, though we discuss an extension of DASh to 
w x = p x /p x ^ -1 models. 

3. EARLY ANISOTROPIES 

The dynamical processes at early times (e.g., acoustic 
oscillations of the baryon-photon fluid, Hydrogen and He- 
lium recombination rates and Silk-damping) are governed 
only by u>b, uj m and fh- Photon density matters as well, 
but this is well-determined from the FIRAS measurement 
of the CMB temperature as T = (2.728 ± 0.004) K (95% 
confidence) (Fixsen et al. 1996). Dark energy parameters 
and the curvature radius are irrelevant since the dark en- 
ergy density at early times was negligible (in most models, 
certainly for a cosmological constant) and the curvature 
radius at last-scattering was much larger than the horizon 
at that time. 

The small number of parameters which are necessary for 
fixing the statistical properties of the CMB at early times 
and on small scales led Tegmark & Zaldarriaga (2000) to 
create a high-€ grid of angular power spectra with grid 
parameters, ns, uJb, ^>d and /;,,. Although SIa and Q,k do 
affect the projection of comoving length scales into angular 
scales, they do so in a particularly simple manner. With 
the grid constructed at fiducial values of = Q*k, = 
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f2 A , they obtain C\ for non-fiducial values of the curvature 
via (Wilson 1982): 

Ci(riK,QA,ub,um, fh) = Ci(n* K ,nx,wb,w m , fh) (6) 

where C; = 1(1 + 1)Cj/(2tt), 

i/i = D z ^(n* K ,n* A )/D z ^(n K ,ci A ) (7) 

and D Z A is the angular diameter distance to z and z pea k is 
the redshift where the visibility function peaks. In section 

5 we derive Eq. 6 and also a version which does not rely 
on any small-angle approximation, as this one does. 

Early anisotropy effects for DASh are also calculated 
via direct numerical solution of the linearized Einstein and 
Boltzmann equations over a grid of parameters. The key 
difference is that DASh stores the Fourier and Legendre- 
transformcd photon temperature perturbation, A;(fc), in- 
stead of Ci . Because of this difference, our grid only needs 
to contain cosmological parameters, and not the primor- 
dial power spectrum parameters. The dimensionality of 
the grid is reduced (and with it the storage requirements) 
and flexibility is increased since we are no longer restricted 
to power-law descriptions of the primordial power spec- 
trum. A typical use of DASh will first take tens of hours 
of computing the A;(fc) grid by a call of CMBfast (Seljak 

6 Zaldarriaga 1996) for each grid point. Only after the en- 
tire grid is computed (we say "pre-computed" ) can DASh 
produce angular power spectra in <~ 1 second, as adver- 
tised. Specifically, the grid is over parameters u>b, u m and 
fh at fixed values of VIk = 1 — fi to t = fi^, = ^a and 
r = 0. For reasons of algorithmic simplicity the current 
implementation requires the number of grid points for each 
grid parameter to be a power of 2. From this grid, we get 
Ci for any u>b, w m5 fh and the primordial power spectrum 
P(k) by performing multi linear interpolation on the grid 
of A;(fc) and then the integral in Eq. 3. DASh can then 
get any Ci, accurate for I > 100, in the model space of {u>b, 
uj d , ^A, P(k)} via Eq. 6. 

Grid boundary and finite grid-spacing effects can be 
minimized by an intelligent choice of the parameters. For 
example, instead of gridding uniformly in uj m we grid uni- 
formly in In Lu rn which makes the interpolation error more 
uniform over the range of Lu m values. The uniformity of 
errors is desirable since if one holds the number of grid 
points fixed a parameterization that has more uniform er- 
rors has a smaller largest error. We discuss variable choice 
more in the next section. 

4. LATE-TIME AND GEOMETRIC EFFECTS 

Although we can use a low-dimensional parameteriza- 
tion of the early anisotropy, the late anisotropy is sensitive 
to additional effects and more cosmological parameters. 
The additional effects are due to geometry, the decay of 
the gravitational potential which occurs in the curvature or 
dark-energy dominated era, and Thomson scattering off of 
the free electrons in the re-ionized intergalactic medium. 

We take two approaches to including these additional 
effects. One approach requires calculation of a second grid 
of A; (fc) (the "low— t grid" which has more dimensions than 
the high-i? A;(fc) grid already described. The other relies 
solely on semi-analytic calculation for the late-time ef- 
fects. The first we will refer to as gDASh and the second as 
sDASh. The sDASh is not completely grid-independent; 



it relies on the high-£ A;(fc) grid as an accurate descrip- 
tion of the photon perturbations at early times and sub- 
curvature scales. Below we first describe the semi-analytic 
calculation of the various effects and then the additional 
grid. 

We model the radiation temperature transfer function 
as being modified by one additive factor and one multi- 
plicative factor: 

Mk) = A] sw (k) + i? ; (r)A* arly (fc). (8) 

where A° arly (fc) is interpolated from the previously de- 
scribed grid, Aj (k) is the late-time contribution from 
the "Integrated Sachs- Wolfe" effect explained below and 
Ri (t) is the reionization damping factor for late-time op- 
tical depth to Thomson scattering, r. The resulting power 
spectrum can thus be written as 

Q = Cf w + Rf (r)C ; carly + Ri ( r )cf w ~ carly (9) 

In the following subsections we describe how we calcu- 
late (7® arly (which gets geometric corrections), Cj sw , the 

c ISW-early term ^ 

4.1. Geometry 

Although the curvature scale is larger than the horizon 
at last scattering, curvature does have effects on the early 
evolution of super- horizon size modes, which are unob- 
servable at the time of last-scattering, but which have 
observational consequences now. That this is the case 
should not be surprising since it is impossible to map, 
without deformation, a space of zero mean curvature onto 
one with non-zero mean curvature. One can see this for- 
mally as a result of the fact that the eigenfunctions of the 
Laplacian are different in spaces of different curvature. 
We will always denote the eigenvalues of the Laplacian 
as k. One can further define a "wavenumber" in curv ed 
space as 1 = k 2 + K/r 2 urv , where r curv = Hq / \/\Qk\ is 
the curvature radius. In the closed case, the spectrum of 
eigenvalues is discrete and /3r curv takes on only integer val- 
ues. Further, /3r curv = 1,2 are pure gauge modes (Bardccn 
1980). The eigenfunctions of the Laplacian in curved space 
are the so-called hyperspherical Bessel functions. We will 
follow the definition and notation of Abbott & Schaefer 
(1986) for the hyperspherical Bessel function and denote 
them as 3>jg(x)- At small distances and short wavelengths, 
$£j(x) = ji(kx) and k = /3. For more details on perturba- 
tion theory and CMB anisotropies in non-flat backgrounds 
we refer the reader to Kamionkowski & Spergel (1994); 
White & Scott (1996). 

We take as our starting point for this calculation the 
Ai(k) already stored in the high-£ grid with fiducial pa- 
rameter values Q\ = f2 A and fix = 0. Recall that this 
grid is for the temperature perturbation today and not on 
the last-scattering surface. Calculating Ci from this grid 
for arbitrary f^A and Qk requires two steps. We must first 
correct for the effects of curvature at last-scattering, and 
then correct for how the projection from last-scattering to 
today has changed. As mentioned, curvature introduces a 
cutoff scale, fc curv , in the spectrum of Laplacian eigenval- 
ues such that fc curv = l/r curv , 0, VS/r curv for K = —1,0,1. 
Our correction for the effect of curvature at the epoch of 
last-scattering is to introduce a cutoff in the integral over 
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fc; i.e., we define 

Cf = (47T) 2 / 



alkk 2 A 9 (k) 2 P{k), 



(10) 



where the g superscript stands for "grid" and implies that 
the quantity in question has been obtained from the grid 
by interpolation. We have found that this simple approxi- 
mation works very well. Note that P{k) here is the power 
spectrum for the flat model. We then use Cf to calculate 
an intermediate angular correlation function C 9 {9) which 
then needs to be stretched to the correct angular diameter 
distance. For monopolc (isotropic) sources emitting from 
a thin shell the shift is particularly simple (Wilson 1982): 

C(9)=C 9 (9'): 

e = 2 £>^(O,O)sin(072), 

sin(0/2) =smh K (Z/2)/D z ^ k (n K ,n A ). (11) 

The function sinhx(x) is defined as sin(x), x, sinh(x) for 
K = 1, 0, —1 respectively. Lcgcndrc transforming the shifted 
C{9) back to £-space then gives us the (7 ; carly f Eq. 9. 

The transformation of the correlation function is only 
exact for monopole sources on a thin shell at fixed red- 
shift. The dipole source due to the peculiar velocities of the 
photon-baryon fluid, the thickness of the last-scattering 
surface, and late-time effects all violate these restrictions. 
We discuss the resulting (very small) errors in section 5. 

At small angles (9 « 1), Eq. 11 can be cast in terms 
of a shift in £, as given by Eq. 6 (see section 5 for a deriva- 
tion). Thus at small angular scales (£ > 10), a simple shift 
in I is sufficient and one need not Legendre transform. 

4.2. Gravitational Potential Decay 
We now show how to calculate the Cj term and the 

^ISW-early crogs tcrm q{ Eq _ g These from ^ 

late-time generation of anisotropy due to gravitational 
potential (metric perturbation) decay. In linear pertur- 
bation theory gravitational potentials are independent of 
time when the Universe is completely flat and matter- 
dominated, but decay in the presence of curvature and/or 
dark energy. As CMB photons pass through the evolv- 
ing potentials, new (secondary) anisotropy is created via 
what is called the Integrated Sachs- Wolfe (ISW) effect (Hu 
& Sugiyama 1995). 

We calculate the ISW effect by evaluating the line-of- 
sight integral: 



Ap w (fc) = f° dn^ws™ 



(12) 



where rjo is the conformal time today and \ — Vo — V- V\&te 
is some late time prior to the onset of curvature domina- 
tion or dark energy domination, which ever is earlier. 
The ISW source term is: 

S ISW (fc,?7) = 2e- T <")*, (13) 

where the optical depth to Thomson scattering is given by 



r(??) 



L 



fair) 



(14) 



where n e is the mean number density of free electrons and 
or is the Thomson cross section. The visibility function is 
g = f exp(— t(t])). The gravitational potential is defined 
in Ma & Bertschinger (1995). 



To calculate "J, when we make the grid of flat models 
we store not only A;(fc) but also ^ 9 (k, z = 100). This we 
use as a transfer function, and then numerically solve for 
the growth factor, D(z) so that 



= ^foL ^ g (fc,100) 



(15) 



D(100) 

This factorization is possible because the evolution of non- 
relativistic matter perturbations is independent of k for 
modes inside or outside the horizon when the clustered 
components are pressureless (Heath 1977). For the growth 
function D(z), we use the approximation given by Carroll 
et al. (1992). 

For values of f3 smaller than some (I dependent) multiple 
of l/n , we evaluate the intergral in Eq. 12 explicitly. For 
other values we use a generalization of the weak-coupling 
approximation of Hu & White (1996) which works best for 
large values of I and (3. Since $^(x) is a rapidly varying 
quantity, one can take the source term out of the integral 
and evaluate it at the conformal time where $^ (x) attains 
its maximum (xmax = Vo — Imax)- This allows Eq. 12 to 
be written as: 



Ap w (fc) 



.-.ISW/ 



(Mmax) / dx&pU), (16) 



P — - 



One is then left with the integral over the hypersphcrical 
Bessel function whose solution can be written as a recur- 
rence relation; we only need the values of the integral for 
I = 0, 1. Denoting the integral in Eq. 16 by 1^, the follow- 
ing recurrence relation can be derived: 

21-1 j|M I - 1 K(l~Ty t1 _ 2 

I - Kl 2 I ' y/p- Kl 2 

(17) 

For open and flat models analytical solutions can be writ- 
ten by taking 770 to infinity; for closed models, we numer- 
ically evaluate 1^ using Eq. 17. 

Now we turn to the ^ ISW_carl y cross term. The largest 
correlation with late ISW comes from the Sachs-Wolfe 
(SW) effect. We currently neglect contributions from the 
primary Doppler and early ISW effects, though including 
them would improve the accuracy. The Sachs- Wolfe radi- 
ation temperature transfer function is given by Ap w (fc) = 
[6 + v I / ] 9 (?7iss) < E > ^(??— »?iss)) where "lss" stands for Last Scat- 
tering Surface taken to be at z — 1100 and [0o + v P] 9 (?7i ss ) 
is the effective photon temperature for the corresponding 
model interpolated from the high-£ grid. With A° arly ~ 
Ap w thus calculated, we then get 



C, 



ISW-oarly 



/■oo 

= (4^) 2 / /3 2 d/3Ap w A° arly P K (/3), (18) 

where Pk{/3) is the curved space initial potential power 
spectrum (Zaldarriaga et al. 1998), and /3k = 0, 0,2/r curv 



for K 



-1,0,1. 



4.3. Reionization 



Thomson scattering smears out our view of the last- 
scattering surface, and therefore damps the early anisotropy. 
This damping is described by the reionization damping 
factor, Ri{t). Our first step to calculating Ri{t) is to ex- 
tract it numerically from models in a grid (pre-computed 
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with multiple calls to CMBfast) over LJb, uj m , ft and r 
with fix = and J1a = 0. These models have no ISW 
effect so we simply set 

R l (Ub,V m , fh,T~) = Cl(u)b,U m , fh 



As we will see later, at high £ Rf 



,T)/Cl(u> b ,UJ m ,fh,0). 

(19) 

= cxp(— 2t), so we ac- 
F ; (l-cxp(-2r)) + 



tually extract Ft instead where Rf 
exp(— 2t). The F; are stored as a function of £/(£ r + 1) 
and we interpolate between the stored values to obtain the 
Fi for a target model with Qk = 0, tt\ = 0. The reion- 
ization multipole is defined by £ r — D*^^ /r) r , where r\ r 
is the visibility function weighted conformal time (Hu & 
White 1997). As pointed out by Hu & White (1997), once 
we have Ri as a function of £/(£ r + 1), it is not changed 
significantly by curvature or A or any other late time ef- 
fect that happens after reionization. The Ri for arbitrary 
model parameters are thus given by Ri(VIk, J^a) = ^f(0, 0) 
where I = l{l r {Q, K , Q A ) + l)/(Z r (0,0) + 1). 

The reionization damping term Rf(r) could also be ob- 
tained semi-analytically. Note that the source terms for 
the early anisotropy (not shown in Eq. 20) all get sup- 
pressed by e _r , independent of £. The ^-dependence of 
the damping factor comes entirely from the non-ISW late- 
time creation of anisotropy at low £ via a source term 

2b k d i / 

9 [aW 



a* 



351 m drj 

where a is the scale factor and b k = (k 2 — 3K/r, 
Thus 



2 

curv 



-2rr<f 



C>(0) 



C{(0) 



(20) 
(21) 



where C ; RI is calculated using S = S m for z < z r ; and 

5 = for z > z r i. 

In writing S Rl we have neglected the Doppler source 
term proportional to the difference in baryon and photon 
fluid velocities. For reasonable values of the baryon den- 
sity, this late-time Doppler effect only starts to become 
important for z ri > 25 (Hu & White 1996). 

We treat reionzation as if it instantaneously occurred at 
z = z T i. We use z r i instead of r, as the input to DASh, since 
z r j is more directly related to observational constraints 
(Becker et al. 2001; Fan et al. 2001). When we sample 
parameter space in order to characterize the accuracy of 
DASh, as described in Section 5, we always keep z r ; < 10. 
For a review of theoretical work on reionization see Loeb 

6 Barkana (2001). 

4.4. A low-£ grid 

We have also implemented in DASh a numerical calcu- 
lation of the gravitational potential decay and geometrical 
effects with the pre-computation of a low £ Af(k) grid. 
This second-grid approach, called gDASh, has the advan- 
tages of speed and tunable accuracy over sDASh. 

The \ow-£ grid, due to its incorporation of the late-time 
effects, necessarily has more dimensions than the high-£ 
grid. We have chosen these extra variables to be f2A/f2 m , 
since this combination controls the ISW effect and ljk 
since this sets the curvature radius. The low-£ grid is less 
sensitive to fh, Ub and ui m than is the case for the high-£ 



grid so we can grid more coarsely in these (Tegmark & 
Zaldarriaga 2000). 

Although Lo m has little effect on late-time generation of 
anisotropy, this parameter directly controls the amount of 
early ISW effect. Since the early ISW effect is not pro- 
jected to us from the last-scattering surface, the angular 
scaling assumed for use of the early grid will introduce 
some errors. Fortunately, these errors are negligible and 
the late grid can indeed be fairly coarse in uj m . 

We join the results of the low-£ and high-^ calculations 
by simply using the low-f calculation up to a limiting value 
/late- Our algorithm for choosing Zi at0 is derived from a 
combination of analytic expectation and experience. First 
we define 



Ik = 404 



(22) 



200 



K 



1/3 



120 



(23) 



with / curv defined for fix > 0, and then we set /i a t e 
max(/ curv , I \). For closed models we set 



max 



25(-^) 1/3 ,U (1-4 



(24) 



which takes into account the effect of both A and curva- 
ture. 2 . 

Our reasoning is that below Zi ate , ISW and curvature 
effects become important. In principal l s = tt/0 s is an- 
other important scale, above which acoustic modifications 
to the intrinsic temperature on the last-scattering surface 
become important. We avoid extending the late grid to 
I > l s because this allows us to grid coarsely in Wb- Fortu- 
nately ii a te is always less than l s . We choose to make the 
switch at /late instead of l s so that he low-£ grid can be 
coarse in ui m ; cutting at higher £ would require finer grids 
in w m to accurately describe the early ISW effect. 

To calculate A;(fc) for the target model from our n- 
dimensional grid we first locate the 2™ grid points of the 
surrounding hypcrcubc. Then for each of these 2" A;(fc) 
vectors we spline-interpolate (and quadratically extrap- 
olate where necessary) on to the k values of the target 
model. This step is necessary because the k values of the 
grid differ from grid point to grid point. A uniform set 
of k values is not desirable since different models have 
different fc-spacing requirements for fixed accuracy spec- 
ification. For closed models a uniform set of k values is 
impossible due to the discrete nature of the spectrum. 

The grid is inefficient (in computing time and storage 
resources) if a lot of the grid points are for models which 
are far from observationally viable. Since we use rectan- 
gular grids, this means we would like to choose parameters 
such that their viable ranges are independent of the values 
of the other grid parameters. A systematic way to do this 
would be to use the eigenvectors of the parameter Fisher 
matrix for some particular experiment (Efstathiou & Bond 
1999). 

We have not pursued this grid efficiency systematically, 
but rather have made the physically motivated choice of 

2 

For closed models in the second-grid approach we we simply set 
'late = max(40, 1a) which works well for the parameter range under 
consideration (|Hk| < 0.3). 
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our low-f grid parameters as u>b, lnu m , ^JVl^/Vl m and u>k- 
Of course, even with this parameterization, we are still free 
to find the parameter eigenvectors. Perhaps doing so will 
further increase the efficiency of the grid and we may in- 
corporate this in future DASh implementations. Not only 
will eigenvectors provide the advantage of a grid with a 
rectangular region of viable models, but there is also an 
advantage in having the well-determined parameter com- 
binations decoupled from the poorly-determined combi- 
nations. Variations in the well-determined ones (over the 
range of their compatibility with data) are most likely to 
reproduce a highly linear response in C; and thus would 
require very few grid points (only two if the response were 
exactly linear). Variations in the poorly-determined pa- 
rameters will generate a non-linear response, but we will 
not need to model these responses as accurately so once 
again will not need many grid points. 

5. ACCURACY 

To characterize the accuracy of DASh we have compared 
a suite of thousands of models calculated using CMB- 
fast and compared them with the same models as calcu- 
lated by DASh. The CMBfast calculations were done with 
very dense fc-spacing, effectively eliminating fc-spacing as 
a source of error. The DASh calculation used the grids 
as described in Table 1. We created our suite of com- 
parison models by first considering all possible models 
with parameter values cut = 0.0145, .018, .022, .028, uid = 
.06, .09, .15, .2, n A = 0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 
n K = 0, ±0.01, ±0.05, ±0.13, ±0.21, ±0.27, z ri = 0,5,10 
and n — 0.95, 1, 1.05 where n is the scalar power-spectrum 
power-law index. Of these models, those within the grid 
boundaries of Table 1 and having 0.3 < h < 1 and a first 
acoustic peak with 160 < I < 280 were placed in the com- 
parison suite. The suite contains 313 models with z r i = 
and Q.k = 0, 937 models with Q,k = and 6,823 models 
total. We refer to the differences between the two cal- 
culating tools simply as 'differences' rather than 'DASh 
errors' because some of the differences are due to errors 
in CMBfast. Also some errors result in no difference be- 
cause they are common to DASh and CMBfast. We have 
not attempted the more challenging task of providing an 
absolute measure of accuracy. 

The results are shown in Figure 2 for sDASh and Figure 
3 for gDASh. In both figures we have shown the maximum 
and average differences for flat models (top panels), for flat 
models with reionization (middle panels), and for flat and 
curved models with reionization (bottom panels). One 
can see that in every case the rms differences are below 
1% for I > 30 for sDASh and I > 10 for gDASh. The rms 
fractional difference (for I < 1000) is in all cases below the 
cosmic variance error of y/TJl. 

The size of the differences relative to cosmic variance is 
encouraging, but unlike cosmic variance the calculational 
errors that give rise to these differences may be correlated 
from t to £. Therefore the target level is closer to l/l 
than y/TJl (see Knox ct al. (1998) and Efstathiou & Bond 
(1999)) since l/l is the cosmic variance error on a band of 
Ci of width /. At all I values, DASh cither meets this more 
conservative criterion or has sub-percent level differences 
with CMBfast. 

DASh errors are certainly small enough to be insignif- 




o.oi L — 1 — 
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Fig. 1. — Differences between sDASh and CMBfast. Solid 
(dashed) lines show rms (maximum) percentage differences. The 
bottom panel is for all 6,823 models in the comparison suite (see 
text). The middle panel is for the 937 flat models in the suite. And 
the top panel is for the 313 flat models with z r i = 0. 
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Fig. 2. — Differences between gDASh and CMBfast. See caption 
for Figure 1. 
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icant for parameter estimation from current data. But 
we do not yet understand the impact of sub percent level 
errors on parameter-determination from all-sky surveys 
such as the Microwave Anisotropy Probe (MAP) 3 . This is 
under investigation. The interpretation of the sub-percent 
level differences with CMBfast is complicated by the fact 
that the CMBfast calculation for some models can have 
percent level errors. 

There are also errors due to neglect of non-linear effects. 
The most important of these are due to gravitational lcns- 
ing and peculiar velocities of the reionized inter-galactic 
medium. We discuss including these effects in the Exten- 
sions section. 

5.1. High I errors 

At high I the top two contributors to error in the DASh 
calculations are interpolation error (at a level of about 
0.5% by grid design) and error from the projection ap- 
proximation. The chief cause of differences in the curved 
case is due to error in our CMBfast calculations arising 
from inaccuracies in the hypersphcrical Bcsscl functions. 

The projection approximation works remarkably well 
and we now turn to understanding that success. This can 
be seen from Figure 3 where we have plotted the differences 
between two flat models, shifted to correct for their 25% 
difference in angular diameter distance to the last scatter- 
ing surface. Both the fiducial and target models in Figure 
3 are flat and were calculated using CMBfast. Even for 
this large a shift, the errors at high I are < 0.5%. We will 
argue that these differences must be due to numerical er- 
rors other than the projection approximation. We expect 
the projection to work just as well for curved models. 

The scaling of 1(1 + 1)C; with angular-diameter dis- 
tance can be derived if one assumes the emission is from 
isotropic sources on an infinitesimally thin last-scattering 
surface. Then the correlation of temperatures at a given 
angular separation is equal to the correlation of temper- 
atures on the last-scattering surface with a given physi- 
cal separation. Therefore one can determine the angular 
correlation function of one model, from the angular cor- 
relation function of another, as long as both models have 
the same physical conditions before and at last-scattering. 
Specifically spatial perturbations to the effective photon 
temperature monopolc, Oq + 5*, with angular correlation 
function C(8) at angular diameter distance r, have cor- 
relation function C (6) — C(0r'/r) at angular diameter 
distance r'. The effect on the power spectrum is to have 
C[ = 1(1 + 1)C;/(2tt) = C v where V = r/r'l since: 

C i = 1(1 + 1) J d(cos6)C(6r' /r)Pi(cos6) 
~l 2 J 6dOC(6r' /r)J (W) 

= I' 2 J dxxC(x)J (l'x) = C v (25) 

where the approximation is accurate for 9 « 1. 

We do not expect the thickness of the last-scattering 
surface to be a significant source of error. We find an- 
alytically that for thickness t and angular diameter dis- 
tance d shifted to d + Sd, the error in C; is of the order 

3 

http: / / map.gsfc.nasa.gov 
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Fig. 3. — Projection errors for different components of the 
anisotropy. The target model is flat with Qa = 0.9, ujf, = 0.02 and 
<^m = 0.f2, while the fiducial model is flat and has U\ = 0. The 
ordinate shows the relative error (in %) for different components of 
the anisotropy produced at or close to the last scattering surface. 
The solid line shows the projection error for the total anisotropy, 
the dashed line for that of the Dopplcr contribution while the dot- 
ted line shows the error in projecting the ISW effect to the target 
model angular diameter distance. 



(6d/d)(t/d) 2 l 2 d 2 Ci/dl 2 which is only important at the sub- 
0.1% level. 

Anisotropies are generated soon after last scattering by 
the early ISW effect due to the decay of the gravitational 
potential in the presence of radiation. Thus for the early 
ISW effect there is a much thicker "last-scattering surface" 
and the projection approximation does not work as well as 
it does for the other early sources. However, the region in 
£ space where the approximation is worst is also where it is 
a highly subdominant contribution to the total anisotropy. 

Note that velocity perturbations on the last-scattering 
surface are a significant contribution to the anisotropy, 
and are not isotropic sources on the last-scattering sur- 
face. Velocity correlations in three-dimensions decompose 
into the correlation between components perpendicular to 
their separation vector and the components parallel to 
their separation vector. At small scales, the radial di- 
rection is nearly perpendicular to the separation vector 
and since the radial component is all that is important 
for the Doppler effect, we are primarily sensitive to the 
perpendicular component of the velocity correlation. This 
component projects like the monopole and the result is 
that dipole sources do not introduce much error, as can 
be seen in Fig. 3. At larger scales the parallel components 
of the velocity correlation also become important, but the 
velocity contribution is sufficiently small at larger angular 
scales that the resulting shifting errors are negligible. 

A comparison of the open models to flat fiducial models 
using CMBfast shows that the high I shifting differences 
are much larger than the errors we expect analytically. 
The 1% to 2% differences are due to errors in the CMBfast 
calculation, presumably in the hypersphcrical Bessel func- 
tions it uses. We compared the CMBfast outputs (at fixed 
oj m and cut) for a flat model, a model with Ok = —0.002 
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Table 1 



Grid Properties 









In w m 


w b 


T 


total 




n k 


ni 


Storage (Mb) 




8 


32 


8 


2 


1 


2,330 


1600 


1600 


22 


570 


High-Z A?(k) 


1 


1 


32 


8 


1 


256 


6000 


4500 


60 


250 


Fi (reionization) 


1 


1 


4 


4 


8 


128 


1600 


1500 


24 


0.1 


range 


0:2.44 


-0.15:0.16 


ln(0.06) : ln(0.28) 


0.01:0.03 


0:0.3 













NOTES. — Entries in the middle rows to the left of the 'total' column show the number of 
values given for each parameter in the various grids. The 'total' column shows the total 
number of models in each grid. The low £ total is less than the product of the other entries in 
the row because models that do not satisfy 0.1 < h < 1.1 and Q, m > and < S1a < 0.9 are 
not calculated. All grids were given the same parameter ranges, shown in the bottom row. 
The High-£ Af (k) and Fi grids were evaluated at fixed Oa = 0.6 and (respectively), wk = 
and t = 0. The last four columns show the maximum k values, number of k values, number 
of £ values calculated and the storage requirements (in Megabytes) for each grid. 



and another with fix = 0.001. (All 3 models had £1a = 0.) 
The differences in the angular diameter distances are less 
than 0.1% and hence we would expect very little differ- 
ence in Ci. Instead we find 1-2% changes to Ci at high 
£ between the flat and curved models, presumably due to 
difficulty in calculating the high I hyperspherical Bessel 
function accurately. Hence we conclude that projecting 
from a flat fiducial model to a curved target model is more 
accurate than the direct calculation using hyperspherical 
Bessel functions. 

Although the current difference plots show we may not 
have met the more conservative criterion of having frac- 
tional errors below 1/Z, we are optimistic that we can get 
there. As data improve DASh will be able to improve along 
with them. Better data mean smaller viable regions of pa- 
rameter space, so the grid boundaries can shrink, allowing 
for greatly decreased interpolation error without increas- 
ing the number of models in the grid. Projection errors, 
already small, will be further reduced since the shifts from 
the fiducial model will be smaller. 

5.2. Low £ errors 

The gDASh can be made arbitrarily accurate at low 
£ by decreasing the parameter grid spacing. This is not 
the case for sDASh where the accuracy is limited by the 
approximations made. Of course, even in the second grid 
method, decreasing the parameter grid spacings requires 
more computational resources, primarily in being able to 
pre-compute and store large number of A/(fc) files. 

The default grid-spacings in DASh have been chosen so 
that the interpolation errors are < 0.5% at high-f . For low 
£, where percent accuracy is not important, the constraints 
imposed on the grid-spacings are more lenient. All the low 
£ sDASh errors are due to interpolation. 

The sDASh errors at £ < 10 are larger than the gDASh 
errors. However, the conservative l/l cosmic variance cri- 
terion is much easier to meet here. Maximum errors of 
the order of 10% can be tolerated. The maximum errors 
always come from the largest values of |Qk|- A large part 
of the error is due to our approximate calculation of the 



cross-term between the early contributions and the late 
ISW contribution. This could be improved as discussed 
earlier. 

Another source of error in the semi-analytic calculations 
at very low Z (Z < 5) is the inability to carry out the projec- 
tion as outlined in Eq. 11. Let r' and r be the coordinate 
distances to the last scattering surface of the target and 
fiducial models respectively. If r' > r, then there is no 
angle in the fiducial model that projects on to 180° in the 
target model. Hence one cannot obtain Ci by a Legendre 
transform. Choosing our fiducial model to be Oa = 0.6 ex- 
acerbates this problem. However that choice for the fidu- 
cial J7a is justified by the aim of getting the high-! spectra 
very accurately for models close to the fiducial one. Note 
that choosing a non-zero f^A to be the fiducial model intro- 
duces another source of error into the low I semi-analytic 
calculations since once needs to subtract off the ISW con- 
tribution due to the fiducial Q\. The error introduced due 
to this is about the same as the error in calculating the flat 
models. One solution to these problems is to implement 
another high-/" grid with Q\ = 0. This increases the pre- 
computing time and the storage requirements. We have 
not implemented this second grid but instead we approx- 
imate the projection effect and tolerate the error in the 
subtraction of the Q\ ISW effect. For the projection, we 
use l'(V + 1) = 1(1 + 1) * r' 2 /r 2 to get to I' = 4 and then 
use the fact that at these large angular scales the spec- 
trum is mostly sourced by the SW effect. Hence the lower 
multipoles can be obtained from I' = 4 by assuming that 
the ratio of C1/C4 is just given by the corresponding ra- 
tio of the integrals over the hyperspherical Bessel function 
squared times the power spectrum. We pre-compute these 
integrals assuming a scale-invariant power spectrum; devi- 
ations from scale-invariance are unlikely to be a significant 
source of error. 

5.3. Speed 

We used a personal computer with Pentium IV proces- 
sors and version 2.96 of the GNU gec g77 compiler for 
all our calculations. Computation of the high-£ A/(fc) 
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grid took ~ 4 hours, of the low-£ A;(fc) grid took <~ 40 
hours and of the Fi grid took ~ 2 hours. After this pre- 
computing, gDASh requires 1.5 seconds and sDASh re- 
quires 2 to 3 seconds to calculate C\. For comparison, 
CMBfast (with the same fc max and rik as used to precom- 
pute our grid — see Table 1) takes 50 seconds for fiat models 
and about 90 seconds for Qk =h models. 

Significant time is spent reading the grid files each time 
a new model is computed. The grids we used (Table 1) 
could be stored in less than a Gbyte of RAM, and doing 
so may speed DASh up by a factor of two. 

6. POLARIZATION AND TENSOR SPECTRA 

The polarization power spectrum is calculated in DASh 
in a manner similar to the temperature spectrum. We 
store the polarization analogs of A;(fc) in the high-€ grid 
and use the same formulae for angular-diameter distance 
shifting and re-ionization suppression. The grid sizes for 
the polarization Aj(fc) are smaller by about a factor of 2 
because polarization A;(fc) is non-negligible over a smaller 
range of k for any given I. The number of grid points 
required to achieve percent level accuracy is comparable 
to the temperature grid. 

Reionization results in a new peak in the power spec- 
trum at very low I (Zaldarriaga 1997). It is the only im- 
portant late-time effect; there is no analogous ISW effect 
for polarization. We have found a simple scaling relation 
that allows it to be calculated rapidly: 

where Cei is auto-correlation of the E (electric or scalar) 
mode of polarization and Cci is the (E mode) polariza- 
tion - temperature cross-correlation. The initial potential 
power spectrum is taken to be P(k) = A (fc//c p i vot ) Tl ~ 4 , 
and Z pivot = (6000 Mpc fe pivDt )/(w m (l + z ri )) 1/2 . The an- 
gular diameter distance shifting is achieved through the 
relation V = l{l' r /l r ) for C E i and V = 1(1' r + 0.5)/ (l r + 0.5) 
for Cci- The fiducial model, denoted by "*", which we 
have used in our scaling relations has r = 0.15, u>b — 0.02, 
u m = 0.16 and f^A = 0.65. Since most of the optical depth 
to reionization is generated prior to the onset of curvature 
or dark energy domination there is no dependence on fifc or 
17a- Also, since the relevant wavelengths are much larger 
than the horizon at last scattering there is no dependence 
on u>b- Eq. 26 is accurate to ~ 10% near the peak of the 
reionization "bump" . When data warrant higher accuracy, 
a low I grid can be created just as in the temperature case. 

For the tensor temperature power spectrum, perturba- 
tions are continually generating temperature anisotropy 
along the line-of-sight so there is no useful late-time/early- 
time split and the simple angular-diameter distance pro- 
jection corrections do not apply. However, a sufficiently 
accurate ( 3%) fit to the tensor temperature power spec- 
trum already exists (Turner 1996) and we have included 
this in DASh. Note that polarization generation from ten- 
sor perturbations only occurs at last-scattering and after 
reionization so these could be included, when data war- 
rant, in a manner analogous to the scalar case. 



7. EXTENSIONS 

There are a number of ways in which DASh could be ex- 
tended. Additional effects can be included such as lensing, 
the Ostriker-Vishniac effect, patchy reionization, gravita- 
tional waves and dark energy with p/p ^ —1. DASh could 
also be extended to calculate the CMB polarization power 
spectra 4 . 

Lensing of CMB photons by mass inhomogeneities on 
our past light cone leads to a smoothing of the power spec- 
trum by a smoothing kernel W\, such that C; = ^2 V W\,Cy . 
This correction can be calculated rapidly and accurately 
as demonstrated by Zaldarriaga et al. (1998) who calcu- 
late it in real space. This smoothing by lensing is a very 
significant effect, leading to corrections at the several per- 
cent level at I = 1000 and tens of percent corrections at 
I = 3000 (Metcalf & Silk 1998). These corrections are 
useful because they break the curvature, Q\ degeneracy 
(Metcalf & Silk 1998). 

Above we have described how reionization generates new 
anisotropy on very large scales as photons pick up some 
of the peculiar momentum of the electrons via scatter- 
ing. Although the linear theory contribution from this 
effect is very small at small angular scales, the second 
order contribution, called the Ostriker-Vishniac (OV) ef- 
fect Ostriker & Vishniac (1986), can be the dominant 
source of anisotropy at I > 3000. The OV contribu- 
tion is Ci ~ 5/iK 2 according to a numerical calculation 
by Springel et al. (2001) and an analytic calculation by 
Ma & Fry (2001); also see recent forecasts for secondary 
anisotropy by Aghanim et al. (2002). Since it is primarily 
a second-order effect it is especially sensitive to the am- 
plitude of the fluctuations and is therefore not sensitive to 
events at high redshift. In particular, since z r j > 6 it is 
insensitive to z r ;. Semi-analytic means of rapidly calcu- 
lating the power spectrum from the OV effect exist and 
could easily be included in DASh (Hu & White 1996; Jaffe 
& Kamionkowski 1998). 

The transition from neutral IGM to ionized IGM is likely 
to go through a "patchy" period in which the two phases 
are spatially mixed. If patches of reionized IGM are small 
enough (comoving linear extent less than about 20 Mpc) 
then the first order contribution to anisotropy from scat- 
tering off of electrons does not experience cancellations 
and can be large (Aghanim et al. 1996; Gruzinov & Hu 
1998; Knox et al. 1998). This contribution from patchy 
reionization is proportional to the redshift width of the 
transition and (l + z r ;) 3 / 2 (Gruzinov & Hu 1998). A small 
z r j means that patchy re-ionization is almost certain to 
be subdominant compared to the non-patchy contribu- 
tion (although see Aghanim et al. (2002) who assume a 
large typical patch size and find the patchy phase may 
contribute significantly) . 

The current implementation of DASh assumes the dark 
energy is a cosmological constant (i.e., p = —p). Whether 
the dark energy is a cosmological constant or not is a mat- 
ter of profound importance in cosmology and fundamen- 
tal physics. Dark energy models based on slowly-rolling 
scalar fields have w = p/p^ —1 and can in principle be ob- 
servationally distinguished from a cosmological constant. 

4 

This has in fact already been done for a preliminary version of DASh 
(B. Gold, private communication) 
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The effect of dark energy on the CMB is solely gravita- 
tional; all we have to do is calculate how the dark energy 
affects <!'. Ma ct al. (1999) have shown that the effect of 
dark energy on the matter evolution can be factorized into 
one function of k and one function of time. We can use 
their fitting formulae for these functions to rapidly calcu- 
late <!' at late times and therefore the ISW contribution to 
anisotropy. 

8. CONCLUSIONS 

We have described a fast and accurate method for calcu- 
lating angular power spectra from the parameters of adi- 
abatic models. Our implementation, DASh, is publicly 
available. Extensions for polarization, tensor modes and 
lensing are straightforward. 

The speed of DASh will be useful for parameter esti- 
mation from CMB power spectrum data which typically 
requires more than hundreds of thousands of models to 
be calculated. DASh is particularly advantageous com- 
pared to grid-based C\ calculations for parameter estima- 
tion techniques which require calculation of C'i at random 
points in the parameter space, such as the MCMC method 
used in Knox ct al. (2001) and described in Christensen 
et al. (2001). A preliminary version of DASh, valid only 
for flat models, was used in Knox et al. (2001). Since we 
pre-compute a transfer function instead of Ci DASh's ad- 
vantage will be greatest for applications with large num- 
bers of primordial power spectrum parameters — such as 
attempts to reconstruct this spectrum from data (Wang 
& Mathews 2000). 

We are grateful to U. Seljak and M. Zaldarriaga for mak- 
ing CMBfast publicly available, S. Colombi, M. Tegmark 
and M. Zaldarriaga for useful conversations, and A. Kosowsky 
for his hyperspherical Bessel function routines. We used 
CAMB (Lewis et al. 2000) many times as we worked through 
various accuracy tests. LK thanks the IAP for their hos- 
pitality during early stages of this work. This work was 
supported by NASA grant NAG5-11098. 
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